perm filename SNEW2.F4[DRW,LCS] blob sn#079381 filedate 1974-12-13 generic text, type T, neo UTF8
00100		SUBROUTINE SS
00200		COMMON X(100),Y(100),N,X1(512),Y1(512),S(100),K
00300		DIMENSION A(2,200),B(2,200),C(2,200),I(200),I1(200),I2(200)
00400	C	FIND MAX AND MIN POINTS AND SET SLOPES
00500		ICOUNT=0
00600	710	DO 100 K=2,N-1
00700		IF(Y(K).GT.Y(K-1).AND.Y(K).GT.Y(K+1)) GO TO 20
00800		IF(Y(K).LT.Y(K-1).AND.Y(K).LT.Y(K+1)) GO TO 20
00900		IF(X(K).GT.X(K-1).AND.X(K).GT.X(K+1)) GO TO 30
01000		IF(X(K).LT.X(K-1).AND.X(K).LT.X(K+1)) GO TO 30
01100		GO TO 100
01200	20	S(K)=0
01300		ICOUNT=ICOUNT+1
01400		I(ICOUNT)=K
01500		GO TO 100
01600	30	S(K)=1001
01700		ICOUNT=ICOUNT+1
01800		I(ICOUNT)=K
01900	100	CONTINUE
02000		ICOUNT=ICOUNT+1
02100		I(ICOUNT)=N
02200	C	FIND POINTS SUCH THAT X(K)=X(K+1) OR Y(K)=Y(K+1)
02300		DO 900 K=2,N-1
02400		IF(Y(K).EQ.Y(K+1)) I2(K)=1
02500		IF(X(K).EQ.X(K+1)) I2(K)=1
02600	900	CONTINUE
02700	C	SET THE SLOPES BETWEEN
02800		IA=1
02900		IX=0
03000	110	IF(IX.EQ.ICOUNT) GO TO 200
03100		IX=IX+1
03200		IA=IA+1
03300	120	IF(IA.EQ.I(IX)) GO TO 110
03400		S(IA)=(Y(IA+1)-Y(IA-1))/(X(IA+1)-X(IA-1))
03500		IA=IA+1
03600		GO TO 120
03700	200	CONTINUE
03800	C	SET BEGIN AND END SLOPES
03900		SZ=1002001.
04000		K=1
04100		IF(S(K)**2.EQ.SZ.OR.S(K+1)**2.EQ.SZ) CALL S101
04200		CALL S102(X(1),Y(1),X(2),Y(2),S(1),S(2))
04210		IF(S(K).EQ.-1001.) S(K)=1001.
04220		IF(S(K+1).EQ.-1001.) S(K+1)=1001.
04300		K=N-1
04400		IF(S(K)**2.EQ.SZ.OR.S(K+1)**2.EQ.SZ) CALL S101
04500		CALL S102(X(N),Y(N),X(N-1),Y(N-1),S(N),S(N-1))
04510		IF(S(K).EQ.-1001.) S(K)=1001.
04520		IF(S(K+1).EQ.-1001.) S(K+1)=1001.
04600	C	RESET THE SLOPES FOR STRAIGHT LINES
04700		DO 610 K=1,N-2
04800		U1=(Y(K+1)-Y(K))/(X(K+1)-X(K))
04900		U2=(Y(K+2)-Y(K+1))/(X(K+2)-X(K+1))
05000		IF(ABS(U1-U2).GT..0001) GO TO 610
05100		I1(K)=1
05200		I1(K+1)=1
05300		S(K)=U1
05400		S(K+1)=U1
05500		S(K+2)=U1
05600	610	CONTINUE
05700	C	ADD POINTS
05800		K1=N-1
05900		DO 300 K=1,K1
06000		IF(I1(K).EQ.1) GO TO 300
06100		IF(I2(K).EQ.1) GO TO 300
06110		FLG=0.
06120		IF(S(K).EQ.0..AND.S(K+1).EQ.0.) FLG=1.
06130		IF(S(K).EQ.1001..AND.S(K+1).EQ.1001.) FLG=2.
06140		IF(FLG.EQ.1..OR.FLG.EQ.2.) GO TO 205
06200		SZ=1002001.
06300		IF(S(K)**2.EQ.SZ.OR.S(K+1)**2.EQ.SZ) CALL S101
06400		U=(Y(K)-Y(K+1)-S(K)*(X(K)-X(K+1)))/(S(K+1)-S(K))
06500		AX=X(K+1)+U
06600		AY=Y(K+1)+S(K+1)*U
06700		XA=X(K)
06800	1610	XB=X(K+1)
06900		IF(XA.LE.XB) GO TO 202
07000		XA=X(K+1)
07100		XB=X(K)
07200	202	YA=Y(K)
07300		YB=Y(K+1)
07400		IF(YA.LE.YB) GO TO 204
07500		YA=Y(K+1)
07600		YB=Y(K)
07700	204	CONTINUE
07800		IF(AX.GE.XA.AND.AX.LE.XB.AND.
07900	     1	AY.GE.YA.AND.AY.LE.YB) GO TO 290
08000	205	K1=K1+1
08100		DO 210 K2=K+1,N
08200		K3=N+K+1-K2
08300		X(K3+1)=X(K3)
08400		Y(K3+1)=Y(K3)
08500		I1(K3+1)=I1(K3)
08600		I2(K3+1)=I2(K3)
08700	210	S(K3+1)=S(K3)
08710		IF(FLG.EQ.1..OR.FLG.EQ.2.) GO TO 280
08800		N=N+1
08900		Z0=X(K)**2-X(K+2)**2
09000		Z1=Y(K)-Y(K+2)-S(K)*(X(K)-X(K+2))-((X(K)+X(K+2))/2.
09100	     1	-X(K))*(S(K)-S(K+2))
09200	     	Z2=X(K)**3-X(K+2)**3-3*(X(K)-X(K+2))*X(K)**2
09300	     1	-1.5*(X(K)+X(K+2))*Z0+3*X(K)*Z0
09400		AZ=Z1/Z2
09500		BZ=(S(K)-S(K+2)-3*(X(K)**2-X(K+2)**2)*AZ)/
09600	     1	(2*(X(K)-X(K+2)))
09700		CZ=S(K)-3*AZ*X(K)**2-2*BZ*X(K)
09800		DZ=Y(K)-AZ*X(K)**3-BZ*X(K)**2-CZ*X(K)
09900		X(K+1)=-BZ/(3.*AZ)
10000		Y(K+1)=AZ*X(K+1)**3+BZ*X(K+1)**2+CZ*X(K+1)+DZ
10100		S(K+1)=3*AZ*X(K+1)**2+2*BZ*X(K+1)+CZ
10200		IF(Y(K+1).LT.YB.AND.Y(K+1).GT.YA) GO TO 290
10300		X(K+1)=(X(K+2)+X(K))/2.
10400		Y(K+1)=(Y(K+2)+Y(K))/2.
10500		S(K+1)=0.
10510		GO TO 290
10520	280	N=N+1
10530		IF(FLG.EQ.2.) GO TO 282
10540		S(K+1)=2.*(Y(K+2)-Y(K))/(X(K+2)-X(K))
10550		GO TO 284
10560	282	S(K+1)=(Y(K+2)-Y(K))/(2.*(X(K+2)-X(K)))
10570	284	X(K+1)=(X(K+2)+X(K))/2.
10580		Y(K+1)=(Y(K+2)+Y(K))/2.
10600	290	IF(S(K).EQ.-1001.) S(K)=1001.
10700		IF(S(K+1).EQ.-1001.) S(K+1)=1001.
10800	300	CONTINUE
10900	C	FIT THE POINTS WITH N-1 CURVES
11000	305	DO 400 K=1,N-1
11100		IF(I1(K).EQ.0) GO TO 310
11200		A(1,K)=0.
11300		A(2,K)=0.
11400		B(1,K)=X(K+1)-X(K)
11500		B(2,K)=Y(K+1)-Y(K)
11600		C(1,K)=X(K)
11700		C(2,K)=Y(K)
11800		GO TO 400
11850	310	IF(S(K)**2.EQ.SZ.OR.S(K+1)**2.EQ.SZ) CALL S101
11900		B(1,K)=2*(Y(K+1)-Y(K)-S(K+1)*(X(K+1)-X(K)))/
12000	     1	(S(K)-S(K+1))
12100		B(2,K)=S(K)*B(1,K)
12200		A(1,K)=X(K+1)-X(K)-B(1,K)
12300		A(2,K)=Y(K+1)-Y(K)-S(K)*B(1,K)
12400		C(1,K)=X(K)
12500		C(2,K)=Y(K)
12600	400	CONTINUE
12700	C	CALCULATE 512 POINTS
12800		M=N-1
12900		R=M/511.
13000		T=-R
13100		DO 500 L=1,511
13200		T=T+R
13300		KT=T+1
13400		KU=T
13500		TT=T-KU
13600		X1(L)=A(1,KT)*TT**2+B(1,KT)*TT+C(1,KT)
13700	500	Y1(L)=A(2,KT)*TT**2+B(2,KT)*TT+C(2,KT)
13800		X1(512)=X(N)
13900		Y1(512)=Y(N)
14000	600	RETURN
14100		END